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We study variable-rate linear quenches in the anisotropic Heisenberg (XXZ) chain, starting at 
the XX point. This is equivalent to swithcing on a nearest neighbour interaction for hard-core 
bosons or an interaction quench for free fermions. The physical observables we investigate are: 
the energy pumped into the system during the quench, the spin-flip correlation function, and the 
bipartite fluctuations of the z component of the spin in a box. We find excellent agreement between 
exact numerics (infinite system time-evolving block decimation, iTEBD) and analytical results from 
bosonization, as a function of the quench time, spatial coordinate and interaction strength. This 
provides a stringent and much-needed test of Luttinger liquid theory in a non-equilibrium situation. 
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While it is difficult to study genuine non-equilibrium 
dynamics in solid state systems due to the presence of 
many relaxation channels (phonons, impurities, interac- 
tions etc.), cold atoms in optical lattices provide an ideal 
laboratory for non-equilibrium investigations due to the 
high degree of control over various dissipation mecha- 
nisms. Cold-atom experiments in the past decade have 
explored a wide variety of non-equilibrium quantum dy- 
namics in previously inaccessible regimes [lj, [2[ . This has 
also led to an increasing amount of theoretical activity 
0, 13 ■ Key issues include thermalization as well as equi- 
libration and their relation to intcgrability [2[ , pumping 
beyond the adiabatic limit or quantum fluctuation rela- 
tions Q , and universal near-adiabatic dynamics in quan- 
tum critical systems 0, Linear quenches occuring 
over a finite time can interpolate between the more fa- 
miliar limits of an instantaneous quench and an adiabatic 
sweep. Very recently, a few experiments have examined 
the response of many-body experiments to such finite- 
time quenches [1, Q . It is thus of vital current interest 
to address the dynamics under linear sweeps of system 
parameters such as interaction strength. 

The response of a system to an external perturbation 
depends sensitively on its spatial dimension, as famously 
demonstrated in the experiment of Ref. Q- There, one 
dimensional interacting bosons did not reach thermaliza- 
tion within the experimental timescale, while their higher 
dimensional realizations did. One dimensional systems 
are notoriously strongly correlated due to the limited 
phase space for scattering. The non-interacting ground 
state is immediately destroyed by interactions, forming 
a Luttinger liquid (LL) in many instances 0, Q , and de- 
scribed by critical phenomena of collective modes with 
anomalous (non-integer) power-law dependence of corre- 
lation functions. 

Quantum quenches in Luttinger liquids have been ad- 
dressed by several authors (icil- 17|. However, it is not 



clear to what extent the Luttinger liquid (LL) picture, 
which is a genuine low energy descri ptio n, is applicable 
under non-equilibrium circumstances [17| . For an abrupt 
interaction ch ang e, certain observables revealed universal 
LL behaviour |15j. While in equilibrium, the relevance or 
irrelevance of a given process can be classified (e.g. using 
power counting), such an approach is not reliable out of 
equilibrium, where additional energy scales emerge (e.g. 
quench duration (l3j and the difference between initial 
and final parameters). 

Therefore, to understand the applicability of the con- 
tinuum LL description for quenches, one needs to go be- 
yond the LL paradigm by either considering additional 
terms in the Hamiltonian (termed irrelevant in equilib- 
rium) or by comparing the results of the LL theory to 
numerical simulations on lattice models. We have under- 
taken the second option, and performed extensive numer- 
ical simulations of arbitrary rate quenches on the XXZ 
Heisenberg model and compared these to bosonization 
results. Similar approach was undertaken in Ref. flBj for 
the case of a sudden interaction quench. 

The model under study is the XXZ Heisenberg model, 
which reads as 
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where m indexes the lattice sites with lattice constant 
set to unity, and J > is the antiferromagnetic exchange 
interaction. We are going to manipulate J z as a function 
of time as J z (t) = J z Q(t), with Q(t) encoding the explicit 
quench protocol, switched on at t = 0. We concentrate 
on a linear quench, namely with Q(t < r) = t/r and 
Q{t > t) = 1. Via a Jordan- Wigncr transformation 
the XXZ Heisenberg Hamiltonian maps onto spinless ID 
fermions with nearest neighbour interaction 
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up to an irrelevant shift of the energy. The c's are 
fcrmionic operators. Alternatively, Sf~ acts as a hard core 
boson creation operator to site I, and the model maps 
to the hopping problem of hard core bosons, interacting 
with nearest-neighbour repulsion. This can conveniently 
be treated using Abelian bosonization, after going to the 
continuum limit as [a, M 
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with u(q) = v\q\ (v = J being the bare "sound velocity"), 
and b q the creation operator of a bosonic density wave, 
and g(q,t) = g2{q)\q\Q{t), g 2 {q) = g 2 exp(-_R |g|/2) with 
Ro the range of the interaction. The connection between 
the two models is established as —1 <C <?2/2u = Jz/nJ "C 
1. The velocity renormalization 0] by J z is neglected 
since it does not affect the physics we discuss to leading 
order. 

We describe time-evolution using the Heiscnbcrg equa- 
tion of motion, leading to 



b q (t)=u q (t) b q (Q)+v* q (t) 61,(0), 



(4) 



where all the time dependence is carried by the time de- 
pendent Bogoliubov coefficients u q (t) and v q (t), and the 
operators on the r.h.s. refer to non-interacting bosons 
before the quench. The bosonic nature of the quasipar- 
ticles requires |w g (i)| 2 — |w g (i)| 2 = 1. The Bogoliubov 
coefficients arc determined from 1131 
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with the initial condition u q (fi) = 1, v q (0) = 0. 

We now obtain various dynamical quantities using the 
bosonization approach and compare the results to nu- 
merical simulation of the quench on the lattice system of 
Eq. |T]). The numerical simulations were performed us- 
ing a combination of a matrix-product state (MPS) [l9[ 
based infinite density matrix renormalization (iDMRG) 
[20l - |22j and the infinite time evolving block decimation 
(iTEBD) [H] algorithms 0, 0] . In our implementation 
of the two algorithms we use infinite, translationally in- 
variant systems. Working in the limit of infinite systems 
has the advantage that no finite size effects show up and 
the only approximation is the finite bond dimension (%) 
of the matrix-product state (MPS). In critical systems 
(as the one we are studying), the finiteness of \ induces 
a finite correlation len gth £ cx x K with k being a model 
specific parameter pBI |26| . For our simulations, we use 
MPS's with bond dimensions of up to \ = 2000 to en- 
sure that the induced correlation length does not affect 
our results. We first use the iDMRG method to find the 
ground state by optimizing variationally a wavefunction 
in the MPS representation. Then the actual quench is 
simulated using the iTEBD technique. This technique 



is based on a Suzuki- Trotter decomposition of the time- 
evolution operator and provides an efficient algorithm to 
perform the real-time evolution of the MPS during the 
quench. We choose a time-step of 5t = 0.01J -1 and use 
a second-order Trotter decomposition. 

The most obvious quantity to start with is the heating, 
i.e., the energy pumped into the system in excess of the 
final ground state energy. This is found to be 
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Here we introduced the microscopic time scale, tq = 
Ro/2v, and E gs = — Lg\j \-kvR\ is the adiabatic ground 
state energy shift to lowest order in 172 , with L the system 
size. The heating in the near-adiabatic limit in ID gap- 
less systems has been addressed in Refs. 27, 28|, where 
non-universal behaviors were reported. The universal 
ln(r)/r 2 heating seen in Eq. ^ was mentioned previ- 
ously in Ref. [29(. In Fig. [TJ we compare Eq. ^ to 
the numerical result, using Rq = 0.5622 as the only free 
parameter, what we obtain from the block fluctuations 

(Fig. ED. 




FIG. 1. (Color online) The heating is plotted from iTEBD 
for J z /J = 0.1 (blue circles), 0.2 (red squares) and 0.4 (green 
triangles) together with the prediction of Eq. © (black solid 
line), using R = 2vt = 0.5622 (from Fig. O J z /J = 0.1 
data). The agreement remains excellent for small variations 
of Ro as well. 

The fluctuations of S z in a given box with size I (in 
units of the original lattice constant) are characterized 
by the quantity 



F(l,t) 
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In the fermionized picture, this corresponds to density 
fluctuations in a box, which was shown to scale iden- 
tically to the entanglement in equilibrium [30j ]. After 
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bosonization, it is given by 

F(l,t) = ^(lcf>(l,t)-Cf>(0,t)} 2 ) = 
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Here, the constant term (1/2) on the r.h.s gives the fluc- 
tuations of free fermions or hard core bosons as F(l) ~ 
ln(7), while the other terms, depending on the Bogoliubov 
coefficients, contain information on the quench. The 
dominant contribution comes from the He[u q (t)v*(t)] 
term to the exponent (l8j . which is of order g^jv, while 
the v q {t) 2 goes only with (^/f) 2 , and can be neglected in 
the perturbative regime. From this, we obtain for t = r 
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where f(y) = \ ~£, s =±i s (v ~ s ) ln \v - s \ and F ( l > °) = 
F Jz=0 (l) ~ hi(Z) accounts for the fluctuations in the ini- 
tial system, which is subtracted to focus on the effect of 
the quench. In the steady state, F(l, 0) — F(l, t — > 00) = 
(<?2/^7r 2 ) hi (l/Ro) becomes independent of the quench 
time, and coincides with the equilibrium fluctuations. 
Right at the end of the quench with t = t ^> l/v, F(l, 0) — 
F(l, t) takes the same value as in the steady state, while 
for vt <C I, it saturates to (32 /W 2 ) ln (r/er ). These an- 
alytical results are compared to the data obtained by the 
iTEBD in Fig. [2J The numerical results for the fluctu- 
ations are more sensitive to truncation effects due to a 
finite matrix dimension x than the other observable we 
calculated. To obtain unbiased data, we performed sim- 
ulations up to x — 2000 and extrapolated the data to 
X^co. For each final J z , the only global fitting param- 
eter for all r's is the short distance cutoff, Ro, which is of 
the order of 0.5 (in units of the original lattice constant). 
For small J z , the agreement is excellent. Remarkably, 
the semi-quantitative agreement persists for J z values as 
high as J z = 0.4. Also, for a given J z , the agreement 
gets better with r, i.e. moving towards the adiabatic 
limit, since the larger r, the smaller the energy pumped 
into the system and the more reliable the Luttinger liquid 
description is. 

The quantities considered so far could have in princi- 
ple been obtained by using adiabatic perturbation theory 
[31| , since our perturbative results capture only the low- 
est order correction in J z /J to the above physical quanti- 
ties. Therefore, we now focus on the spin flip correlation 
function (S + S~), which contains the bosonic fields in the 
exponent and demonstrates the non-perturbative nature 
of bosonization: the present approach yields to first cor- 
rection in J z to the exponent of the spin-flip correlation 
function. Perturbation theory would only yield the low- 
est order correction to the whole correlator, and not to 
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FIG. 2. (Color online) The fluctuation of S z in a box of 
length I is plotted for J~/J = 0.1 (left panel), 0.2 (center 
panel) and 0.4 (right panel) for Jr = 3.8, 8.7, 20.2 and 71.2 
from bottom to top. Points are iTEBD results and lines are 
fits using the bosonization result Eq. (|10[) . The only global 
fitting parameter in each panel is determined as Ro = 0.5622 
(left), 0.6275 (middle) and 0.7446 (right) in units of the lattice 
constant. 



its exponent. Thus, it requires the non-perturbativeness 
of bosonization to account for the numerical data and to 
produce power-law correlation functions. 

The most singular, staggered part of the transverse 
magnetization [9j is given by 
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(11) 



where 9{x) is similar to Eq. ((5J), except for an extra 
sgn(g) multiplier within the q summation (l8j . S + (x) is 
also the hard core boson creation operator in the contin- 
uum limit, and a is a short distance cutoff. 

The spin flip correlation function of the XXZ model, 
which corresponds to the hard-core boson single particle 
density matrix, Gs{x,t) = (S + (x,t)S~ (0,i)) is obtained 
as 
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where {[0(x,t) - 8(0, t)} j is similar to Eq. ©, only the 

sign of the last term is flipped [HI . 

Right after the quench at t = r and in the |x|,«r 3> Ro 
limit, the spin flip correlation function reads as 
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where C = \fe2 6 stems from the correlator of 

hard core bosons on a lattice in, e.g., the XY model (52 = 
0), A = 1.28243... is Glaisher's constant [32|. These 
non-perturbative results are tested against numerics in 
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Fig. |3l where we fix R = 0.5622 from Fig. [2j Similarly 
to the previous comparisons, the agreement is excellent 
and works qualitatively upto rather large J z . 




FIG. 3. (Color online) The spin flip correlation function is 
shown for J z /J = 0.1 (left panel), 0.2 (middle panel) and 0.4 
(right panel) for Jt = 0, 20.2 and 108.3 from top to bottom 
with Ro = 0.5622 (from Fig. [2j from Eq. (O, together with 
the numerical data. The power-law exponent changes from 
— | — ^fj for x <C vt to — | for x ^> vt, as \G(B(x, r)\y /f \x\ « 
C(R /mm[x, 2«r/e]) 92/2u . Results from the XY model fix 
the prefactor of the correlation function as well, leaving Ro as 
the only adjustable parameter. The r = results correspond 
to that in the XX Heisenberg model [32]. At short distances, 
the correlator is strongly influenced by the presence of the 
lattice. 

The short distance behaviour (< vt) in Figs. [5] and 
|3] is dominated by high energy (> 1/t) modes, evolving 
adiabatically. The correlators thus behave identically to 
the adiabatic case (r — > oo). However, the long dis- 
tance (> vt) response is dictated by low energy (< 1/t) 
modes, feeling a sudden quench, and the observables in 
this range reveal the sudden quench behaviour (r — > 0). 
We have also checked that the numerical data for time 
dependent correlators are also successfully described by 
our bosonization scheme. 

After the quench (t > r), Eq. (|12|) still applies after 
changing r to t. The momentum distribution (MD), i.e. 
the spatial Fourier transform of Eq. (fT2j) , to first order 
in g 2 behaves as 
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where k — \ \k\ — ir\. In the steady state, it remains identi- 
cal to the adiabatic expression [8|,|9( in spite of the quench. 
Had we taken a ferromagnetic coupling (J < 0), the di- 
vergence would occur at k = as is the case normally 
for hard core bosons j33|. The steady state (t oo) 
response thus coincides with the equilibrium one to first 
order in the exponent, irrespective of the quench time. 
Hig her order terms, however, will modify the exponent 
13| . Eq. (fT4| is directly accessible experimentally using 



librium purview of this description, by deriving quan- 
tities using an out-of-equilibrium Luttingcr liquid the- 
ory and comparing them to exact numerical calculations 
using iDMRG/iTEBD for the XXZ chain. Since sev- 
eral calculations have appeared in the literature treating 
the Luttinger model in non-equilibrium situations, it is 
important to develop intuition for the reliability of the 
Luttingcr model as a description of the non-equilibrium 
physics of lattice models. Our work is an important 
step in that direction (cf. Ref. [Hj]). Remarkably, even 
though our bosonization calculations arc pcrturbative in 
J z , they provide an excellent quantitative description 
even for moderately large J z values. 

Our work opens up a number of new questions worth 
pursuing in future research. We have found bosoniza- 
tion to describe well linear-quench dynamics from J z = 
upto moderate values of J z . While this is indicative of 
the broad applicability of bosonization out of equilib- 
rium when starting from an initial ground state, it might 
also be fruitful to explore similar issues for other non- 
equilibrium situations. In particular, one might wonder if 
the Luttinger model is quantitatively useful for instanta- 
neous quenches involving large changes in J z beyond the 
observables considered in Ref. , or for cases where the 
initial state is not a ground state. In general, it is not 
well-understood which non-equilibrium situation might 
make which type of irrelevant or marginal operators im- 
portant. 
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time-of-flight imaging of quenched hard core bosons. 

To summarize, we have applied the Luttinger model 
description for a lattice model outside the usual equi- 
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EPAPS supplementary material 

Here we detail and explain some technical aspects of the calculation, presented in the main text. The Jordan- Wigner 
transformation reads as [8f 
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where the c's are fermionic operators. 

In terms of the time dependent Bogoliubov coefficients, the heating reads as 
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where both terms contribute in (g2/v) 2 order. 
For a linear quench and t > r, 



v(q,t) = [e>q>(»w(g)(t - 2r)) - exp(iu(q)t) + 2iw( 3 )re3q)(-«j(g)t)] 
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and 



u(q,t) = exp(— iuj(q)t) 
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to lowest order in g 2 . Consequently, 



Rc[u q (t)v* q (t)} 



.92(g) 
2v 



sm(2v\q\(t - r)) - sin(2v\q\t) 
2v\q\r 



(A.19) 



to lowest order in gijv. In the steady state (t — > 00), the trigonometric functions average to zero and only the 
—g 2 (q)/2v factor determines the response. The adiabatic limit, r — > 00 leads to identical result, therefore the steady 
state behaves as the adiabatic limit, as if no quench occured to the system. 
The bosonic field, 9{x) is defined as 



) = E^ g %)(^- Qkl/ % + ^-) 



(A.20) 
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ans its autocorrelation function reads as 

([6(x, t) - 9(0, t)} 2 ) = E sin2 (f ) (\ + KWI 2 - Mu q (t)v* q (t)]) . (A.21) 

A quantum state on a chain of length L can be written in the following MPS form: 

|*)= A^A^...A^\j u ...,j L ), (A.22) 

h,—,3L 

where are \n x Xn+i matrices, and \j n ) represent local states at site n. The matrices at the boundary (i.e., 

n = 1 and L) are vectors because the outer index is zero dimensional. The MPS representation is efficient in one 
dimensional systems because it exploits the fact that the ground-state wave functions are only slightly entangled [IH ■ 
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